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The  purpose  of  this  thesis  is  to  develop  an  algorithm  Dr  segmenting  images 
corrupted  by  a  high  level  of  noise  with  different  characteristics.  In  particular  the 
images  considered  are  composed  of  several  regions  describing  different  objects  and 
background.  The  algorithm  described  is  based  on  a  Markov  Random  Tie!  J  t MRl  : 
model  of  the  image  and  uses  Kalman  filtering  (Kid  techniques  and  Dynamic 
Programming  (DP)  in  order  to  smooth  within  the  regions.  The  theoretical  background 
for  one  dimensional  and  two  dimensional  data  which  have  dilferent  characteristics  and 
simulation  results  are  presented,  with  examples  of  synthetic  data  and  underwater 
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I.  INTRODL  CTION 

The  objective  of  segmentation  is  to  divide  a  given  image  into  meuninglul  regions 
or  units.  In  many  vision  problems  the  first  task  is  to  distinguish  objects  to  discrnunate 
between  various  obstacles. 

A  general  vision  system  for  Artificial  Intelligence  applications  can  be  decomposed 
in  the  following  blocks: 

1.  Data  acquisition:  The  picture  of  interest  can  be  taken  by  suitable  cameras  lor 
these  purposes  and  can  be  digitized  using  proper  software  and  hardware. 

2.  Segmentation:  Removing  noise  and  dividing  the  image  into  significant  regions. 

a.  Image  understanding:  Deciding  which  objects  are  present  and  classifying  them 
according  to  size,  shape,  etc. 

The  general  approach  at  the  basis  of  image  segmentation  is  based  on  the 
identification  of  different  properties  characterizing  the  regions  of  interest.  Tor 
example,  several  objects  can  be  identified  by  their  average  intensity  levels  or  by  the 
textures  on  the  surfaces.  The  task  of  an  image  segmentation  algorithm  is  therefore  to 
identify  the  regions  from  the  vision  signals  which  contain  noise  (due  to  the  electronic 
equipment  and  other  environmental  factors  such  as  turbid  waters  in  underwater 
environments)  or  uninteresting  details.  For  example,  if  it  is  desired  to  recognize  the 
presence  of  a  house  in  the  picture,  we  need  to  segment  the  image  by  disregarding  the 
doors,  windows,  cracks  on  the  wall.  etc. 

Although  all  the  methods  in  the  literature  for  segmentation  are  well  suited  to 
most  of  the  cases  of  interest,  they  have  poor  performances  in  the  presence  of  a  high 
level  of  noise.  In  the  cases  of  noisy  data  statistical  techniques  are  moie  suitable.  In 
particular  these  are  based  on  classical  techniques  of  estimation,  such  as  Maximum 
Likelihood  (ML)  or  Maximum  a  Posteriori  (MAP).  Statistical  estimation  techniques 
rely  on  the  use  of  distinct  statistical  models  for  the  original  image  and  lor  the 
disturbances. 

In  this  thesis  we  present  a  segmentation  algorithm  for  images  affected  by  a  high 
level  of  noise,  based  on  statistical  techniques.  In  particular  we  use  Markov  Random 
Fields  (MRF)  as  a  model  of  the  original  image,  and  we  assume  White  Gaussian  noise 
as  a  model  for  the  disturbance.  The  leason  for  the  choice  of  MRF  is  based  on  several 
considerations: 

1.  They  extend  one  dimensional  well-known  models  to  the  multidimensional  case. 
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2.  They  relate  global  statistics  to  local  dependencies  of  the  data. 

3.  They  can  model  compact,  distinct  regions  by  a  suitable  choice  of  parameter';. 

The  smoothing  is  obtained  by  a  combination  of  MR!'  (for  the  transition  between 

regions i.  Kalman  Filtering  Techniques  (to  filter  within  the  regions)  and  Dynamic 
Programming  (to  determine  the  best  sequence  of  edges  which  maximizes  the  likelihood 
function ). 

Related  to  the  works  of  segmentation,  the  algorithms  of  Derin.  Flliott  and  Cristi 
[Refs.  1.2).  Geman  [Ref.  3|.  Besag  [Ref.  -!|  base  shown  the  effectiveness  of  these 
techniques  segmenting  images  corrupted  by  a  high  level  of  noise.  Parallel  to  these 
works,  a  combination  of  Autoregressive  and  MRF  models  has  been  presented  by 
Therrien  [Ref.  5j  to  segment  textured  images  of  terrain  data. 

In  the  next  chapter,  the  properties  of  modeling  the  original  scene  by  Markov 
Random  f  ield  and  choice  of  the  parameters  in  Markov  Random  Fields  arc  presented. 
Fstimation  algorithms  for  one  and  two  dimensional  data  are  given  in  Chapter  III. 
Finally,  simulation  results  for  different  characteristics  of  data  are  the  subjects  of 
Chapter  IV. 
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II.  STATEMENT  OF  THE  PROBLEM  AND  MODELING  l  SING  MKi 

I  he  problem  to  be  addressed  in  tins  chapter  is  how  to  model  the  or;..:..,!  .:.c 

and  its  estimation  from  the  given  (observed)  data.  As  stated  in  the  l;.t,o  imt.  the 
model  for  the  original  scene  is  assumed  to  be  a  Markov  Random  Livid  (MRl  .v.,o-c 
properties  are  discussed  in  the  next  section.  As  explained  below,  the  cMimat:  m  o 
based  on  a  Bayesian  approach  and  a  Maximum  a  Posteriori  (  MAI’>  teshmauec 

A.  PROBLEM  STATEMENT 

Let  the  original  scene  X  =  !Xk  :  be  a  random  lit  1J  on  a  Unite  two  dimensional 
lattice  L  =  nk.t)  e  Z;  :  1  <  k  <  N,.  1  <  t  <  N,)  where  Z  is  the  set  of  integers. 

assuming  discrete  values  F,.  F,.  F, . which  are  constant  in  regions  R.,  R„  R  ....  with 

R,  c  L  and  Rj  H  L=  0  for  i*j.  For  example  this  might  be  the  case  ot  sex er  d  cbje-ts 
with  dillerent  intensity  levels. 

Also  let  X  be  corrupted  by  an  additive  noise  and  modeled  by  a  random  held 
'V  =  i\Vkth  (k.t)  e  L.  Furthermore.  \V  will  be  assumed  to  be  a  White  Gaussnm  Field 
with  zero  mean  and  known  variance  <T~.  W  is  identically  independent!}  distributed 
ua.d.i  and  independent  of  X. 

1'heretore.  the  observed  unaue  can  be  eiven  bv  a  random  field  A'  =  :  Y.  ,  as 


Vka  =  gk.t'Xk.f  Wk.t» 


where  (k.t)  e  L  and  gk  t  (,  .  .)  is  an  arbitrary  function.  In  our  case  we  assume 
noise,  and  therefore  the  observed  signal  can  be  expressed  as 


Vt  =  xk.t  +  wk., 


Now,  the  problem  is  to  estimate  X  from  the  data  Y.  Since  we  assume  tli.it  only 
noisy  observations  are  available,  the  estimate  can  be  obtained  by  mnximi/nrioa  of  the 
a  posteriori  distribution  of  the  Scene  with  respect  to  x  by  a  Bayesian  approach. 
Therefore 

rx.v<xly>  =  n\a\,.i  Py:\-ix;v)  .  2.M 


lx:  1  XiV 


» 


m 


where  x  is  the  estimate  of  original  scene  anJ  Py  y<\iy »  =  IVobul\;:r.  .  X  = 

As  it  can  be  seen  from  Equation  2.  \  the  probiem  is  to  niaxina/o  '!.e  .  ; 

,  A  .  ... 

probability  ol  \  given  y.  Inis  form  o!  Bayesian  ^.siimatKn  is  kr;nv. i.  is  \,V 

P.iK.’rnoii  or  MAP  estimation. 

Bayes  factorization  and  the  convenience  u!  logunthmic  u  penmen  yields 

InPy.^ty  \t  -  InPyix)  =  max,^,  :InPy  yiy  \i  -  InPy.  i, 

1  he  two  likelihood  terms  on  the  let t  and  riel:;  hand  sides  of  l.qua:;  . 
determined  by  the  model  of'  X  and  of  tiie  disturbance,  li.e  model  lor  the 
assumed  to  be  a  Markov  Random  I  ield  i.MRI  »  which  is  formally  d.ik.e.i 
| Ref.  3:  p.  “24|  and  by  L'lliot.  Denn,  Crtsti  and  (jeman  [Refs.  <>.“]. 

B.  MARKOV  RANDOM  FIELDS  (MRF) 

Definition:  Let  L  be  a  finite  lattice  L  =  p.k.t)  :  1  <  k  S  \(.  1  ^  t  i 
qk  (  defined  as  a  subset  of  L  so  that  (k.t)  i  qk  (  where  ik.t)  e  L  and  u.ji  e  i|  . 
only  if  (k.tfe  ■.  Then  a  random  field  X  =  |Xk[;  with  the  property 

P:  x,  ,  =  ;  X- .  =  \  . .  i  i.j  i  e  l  ;  = 

k.t  k.t  I.J  1. 1  ' 

P:.\\  =  x,  ,  j  X  .  =  .v  .  .  (i.j >  e  ii.  J 

k.t  k.t  1  I.J  I.J  '  ‘  'K.t' 

is  a  Markov  Random  Field  with  neighborhood  qk  . 

A  Markov  Random  Field  can  be  illustrated  as  in  Figuie  2.1.  where  the  s 
of  the  clement  (k.t)  indicated  by  a  depends  on  its  neighbors  indicated  by 
only.  I  wo  simple  neighborhoods  are  shown  in  Figure  2.2.  These  are: 

Hi  ^  =  1(  l.n)  :  0  <  (i-1  )2  +  ( j-n):  ^  11 

and 

Huj;  =  {(l.n; :  0  <  (i-1)2  -  (j-nr  X  2| 

Relatively  simple  neighborhood  systems  as  in  Figure  2.2  are  adequate  in  modeai 
scenes  of  interest. 
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Figure  2.1  Dependence  of  Markov  Random  Field 
P^i  F Everything  else). 
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(a)  (b> 


Figure  2.2  Neighborhoods  n  1  (ai  and  n.  ‘  (b). 

i.j  '  i,j  ' 


Although  a  MRF  can  be  considered  as  a  multidimensional  extension  of  a 
Markov  Cham,  a  difficulty  exists  in  defining  MRF.  In  fact  the  mam  concept  at  the 
basis  of  Markov  Chains  is  the  concept  of  transition  probabilities.  These  are  subje.  '<  to 


mild  conditions  and  they  rely  on  the  fact  that  an  oidering  can  be  dctert.aned  in  one 
dimensional  problems.  This  is  not  the  case  in  multidimensional  problem'-  w hue  an 
ordering  of  points  in  general  does  not  exist,  which  presents  a  difficulty  in  the  definition 
of  transition  probabilities.  Because  of  this  we  have  to  approach  MRTs  m  a  duicient 
way  from  the  Markov  Chains.  Iherefore.  to  determine  the  '-tructure  for  tire  joint 
probability  of  MRF  models,  we  need  some  new  delinittons  and  theorems.  In 
particular  the  joint  probability  is  based  on  the  concept  of  i-iUnte  gr.en  below: 

Definition:  Given  a  Markov  Random  Field  iMRl  i  with  neighbors  q  a  n.-me  is  a 
set  of  pixels  which  are  neighbors  of  each  other  [Ref.  (>:  p.l9S|.  1  he  various  types  of 
cliques  for  rj,  j*  and  r||  .*  are  shown  in  Figure  2.3.  Based  on  the  definition  of  cliques, 
the  joint  probability  of  a  Markov  Random  Field  with  neighborhood  q  can  be  expressed 
by  the  following  theorem: 

Theorem  (Hammersley  and  Clifford)  [Ref.  6:  pp.  192-199].  Suppose  X  is  an  MRF 
such  that  P\(M  '■  0  for  all  x,  with  neighborhood  r).  Then  Pv(x>  can  be  defned  as 
following: 


px(x)= Ye‘L(x) 


where 


C(x,  =  VVc(x, 


t.  2.  S  • 


i  2.9) 


«  is  the  set  of  all  cliques  as  shown  in  f  igure  2.3  and  usually  L’(.\)  is  defined  as  the 
energy  function.  V'c(x)  the  potential  associated  with  clique  c.  and  the  partition  function 
Z  is  defined  as 

Z  =  Ve-G(M  (2.10, 

x 


which  is  a  normalizing  constant  that  causes  to  be  a  consistent  probability  measure. 
On  the  basis  of  the  above  theorem  we  can  arbitrarily  assign  a  joint  probability  of 
MRFs  by  the  potential  functions  V  (x).  The  only  constraint  on  V  j x )  is  that  it  has  to 
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Definition :  A  MRP  which  has  ail  nonzero  probabilities  as  m  Equation  2. !  1. 
namely,  P^(x)  >  0,  is  called  a  Gibbs  Field  (GP).  The  Gibbs  model  .an  be  iu;d  to 
model  the  spatial  interactions  between  regions,  in  which  a  pixel  and  its  neighbors  have 
similar  statistical  properties.  In  other  words,  it  models  spatial  contimm;.  v.  h:..n 
characterizes  the  images  we  want  to  segment.  Note  that  the  structuie  c!‘  P  v  (  \  i  is 

^  .  v 

relatively  simple  for  t|j  as  in  liquation  2.1 1.  Its  complexity  increases  considerably  ;er 
laiger  neighborhoods.  Por  that  reason,  in  general  onlv  the  lower  sets  ii  . 1  and  n  are 

'  “  ’l.|  *1 -1 

considered  practically  applicable,  and  even  in  ~  cliques  with  only  one  or  t.vo 
elements  are  taken  into  account.  In  addition,  the  Gibbs  distribution  might  a, so  be 
used  to  model  textures.  In  our  problem  we  just  want  to  model  the  tact  that  the  most 
likely  scenes  have  regions  which  arc  clusters  of  pixels,  and  the  Markov  Random  Picid 
model  is  adequate  for  this  purpose,  provided  by  properly  assigning  the  potential 
functions.  A  particular  model  we  are  going  to  consider  in  this  thesis  is  the  MRP  with 
neighbors  r| 1  as  in  Figure  2.3  (a)  and  with  potential  functions  defined  as: 


vi<xk.t>  =  0 


Vxk.t-rxk.t)  =  vi(xk-i.t’xk.t> 


Pifxu  =  Xk.t- 
-P  otherwise 


(2.13) 


for  any  value  of  k  and  t  and  p  a  positive  constant.  The  larger  the  parameter  p.  the 
more  the  joint  distribution  Px  is  peaked  around  smooth  realizations.  By  this  del'nition 
we  can  write  the  joint  probabilities  as 


lnPx(  x)=  P(Vg(  xkt,xk<t.1 )  +  g(  xk  it.xk. u  )i  -InZ 

k.t 


(2.14) 


where 


e(x.  x  )=<rilfXk-c  Xk’!-1 
k.t  '  k-l,t  ^  otherwise 


(2.1.5) 


Also,  by  the  assumption  of  Gaussian  noise  with  zero  mean  and  standard  ue\ union  <r. 
Equation  2.14  with  Equation  2.4  yields  the  likelihood  as 


sV/ 


V 

;?• 

w 

& 


1 


I 

I 

8 


fi 

$ 

% 


$ 

*K 

•2 

iC 

••5 

‘V 


•*»j 


1 


1 


lnPV!X(y|x)  +  lnPx(x)  =  -_V(  ykifxkitr  +  P  ilg(xk>t.xk.,<t>  -  gl  x^.x^  )j  -InZ 


<  2.1oi 


Hence,  given  the  data  y=  fyk  }  in  noise,  the  noise  model  W(0.  cn  and  the  Gibbs  field 
model  Equation  2.14  with  parameter  p.  the  estimate  of  the  original  scene  is  given  by 

A 

x  =  [xk  1  such  that 


^I(>'ky\a):'P‘Is(xk.rxk.t-i)  +  8(xk.rxk-1.t»} 


(2.17) 


is  minimal. 

In  particular  there  is  no  need  to  undertake  the  difficulty  of  calculating  the 
partition  function  Z  since  it  is  a  normalizing  constant.  So.  it  was  omitted  in  Equation 
2.1~.  Approaches  to  minimization  of  Equation  2.17  can  be  based  on  relaxation 
[Ref.  3:  pp.  730-732]  or  deterministic  [Ref.  4:  pp.  266-267]  or  Dynamic  Programming 
[Ref.  2:  pp.  44-45  and  Ref.  S:  pp.  12-13].  In  the  next  chapter  a  different  approach  to 
minimization  of  Equation  2.17  based  on  Kalman  Filtering  techniques  will  be  presented. 

C.  CHOICE  OF  THE  PARAMETER  IN  THE  \1RF 

The  difficulty  in  using  the  Gibbs  Distribution  as  a  region  model  is  to  estimate  the 
parameters  of  the  model  from  specific  realizations.  Some  methods  have  been  used 
previously  to  estimate  the  model  parameters.  For  example.  Besag  [Ref.  6:  pp.  21 1-212] 
suggested  the  coding  method  where  the  parameters  are  determined  from  subsets  of  data. 
This  requires  solution  of  a  set  of  nonlinear  equations.  Another  example  of  the 
parameter  estimation  method  is  given  by  Derin  and  Elliot  [Ref.  2:  pp.  43-44]  which 
uses  standard  linear,  least-squares  estimation.  This  has  been  proved  to  be  effective  in 
modeling  textures  by  GFs.  A  different  approach  was  used  by  Geman  [Ref.  3:  pp. 
727-729]  which  defined  a  simulated  annihiling,  where  stochastic  relaxation  was 
combined  with  monotonic  increasing  parameters. 

For  this  thesis,  the  parameter  P  of  the  model  in  Equation  2.17  was  set  by  trial 
and  error  until  a  reasonable  filtering  was  reached.  The  optimum  value  of  the 
parameter  P  is  smaller  for  high  values  of  signal  to  noise  ratio.  So,  it  is  necessary  to 
modify  this  parameter  as  the  signal  to  noise  ratio  changes. 
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HI.  SEGMENTATION  AND  SMOO  THING  <>|  ONE  DIMENSION  AL  AND 

TWO  DIMENSION  U  SIGNALS 


The  content  of  this  chapter  is  the  segmentation  alguntnms  lor  one  dimensional 
and  two  dimensional  signals.  Following  a  Miudmum  »•  P-jUtri^n  approach,  estimate  of 
the  original  data  is  obtained  by  minimization  ol  the  cost  Inaction  given  by  Equation 
2.1'.  To  provide  this,  a  smoothing  technique  is  dowsed  which  is  based  on  a 
combination  of  Kalman  Filtering  (KF)  and  Dynamic  Programming  (DPi.  Dynamic 
Programming  is  used  to  determine  the  best  sequence  of  edges  which  maximizes  the 
likelihood  function  and  a  detailed  algorithm  is  presented  in  Appendix  A.  Kalman 
Filteiing  is  used  to  filter  within  the  regions.  This  filtering  technique  is  preferred 
because  the  Kalman  Filtering  is  the  best  linear  optimal  estimator.  If  the  noise  m  the 
data  is  assumed  to  be  Gaussian,  the  Kalman  Filter  gives  the  minimum  variance 
estimate  of  the  original  scene.  In  particular  it  evaluates  the  conditional  nv.\m  of  the 
original  data  given  the  past  observations  (measurements).  If  the  assumption  of 
Gaussian  noise  is  removed,  the  KF  yields  the  best  linear  estimate. 

In  the  next  subsection,  the  segmentation  algorithm  for  a  binary  sequence  (i.e.. 
one  which  assumes  only  two  levels)  is  presented  and  the  result  will  be  extended  to 
general  multilevel  eases.  Related  computer  programs  are  included  in  Appendix  B. 

A.  BINARY  SEGMENTATION  USING  DYNAMIC  PROGRAMMING 

ALGORITHM 

Suppose  that  it  is  desired  to  segment  a  binary  sequence  having  logic  levels 
and  '<)"  which  are  related  to  two  intensity  levels  (gray  levels)  corrupted  by  an  additive 
noise.  The  noise  is  considered  as  a  zero  mean  Gaussian  noise  with  a  known  standard 
deviation  a.  In  the  general  image  processing  extension,  we  can  consider  a  logic  "0"  as 
a  background  and  a  logic  "1"  as  an  object.  An  example  of  the  original  sequence  and 
the  noisy  sequence  is  given  in  the  next  chapter. 

Let  x.  be  the  logic  values  of  the  original  data  of  intensity  levels,  which  assigns 
real  numbers  F(0).  F(l)  to  the  regions  denoted  by  logical  zeroes  and  ones  respectively. 
Then  the  noisy  data  can  be  given  as 

Y.  =  F(x.)  +  W.  (ml) 

I  V  l  l 

where  \V.  ~  N  (0  ,  <7 ).  We  can  define  the  signal  to  noise  ratio  (SNR)  as 
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SNR  = 


1 1 1 )  -  r  < <>> 


In  tins  section  we  assume  to  know  the  levels  f(<U.  hi)  and  the  noise  vj;  ianc-c  cr. 
Here,  the  assumption  of  F  and  c  known  is  not  restrictive,  since  various  algorithms 
exist  in  the  literature  which  enable  one  to  estimate  means  and  variance  of  nixt  tires  of 
Gaussian  populations.  We  have  seen  in  Chapter  d  that  a  MRF  can  be  arbitral ily 
assigned  by  the  potential  function  V.  In  particular  we  can  model  spatial  contmuitv  by 
assigning  hitth  probabilities  to  smooth  signals.  I  his  can  be  achieved  by  assigning  the 
likelihood  lunotion  in  the  following  form: 


**v*i . xN,i=ni?(si+l.v— Is;- 


where 


gix.  .  x,)  = 


1  if  x,  =  x, 

-1  otherwise 


The  parameter  p  models  the  smoothness  of  the  original  data.  In  recursive  form 
the  likelihood  function  can  be  written  as: 


V+-i^xu,xi . xk+i)  Mxo-xi . -V  +  Pgtxk  +  rXk)-v  '}‘k+ i Ji* 

with  t_,  =  0  .  The  estimate  x  is  obtained  by  maximizing  C  over  all  possible  realizations 
xr  Applying  Dynamic  Programming!  DPt  algorithm  in  Appendix  A.  the  maximum  of 
fk  is  obtained.  Results  of  segmentation  of  noisy  binary  image  is  given  in  Chapter  4. 
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B.  SMOOTHING  OF  PIECEWISE  CONSTANT  SIGNALS  IN  ONE 

DIMENSION  USING  KALMAN  FILTERING 

The  algorithm  presented  in  the  previous  section  lor  binary  data  is  suitah.e  tor  a 
small  number  of' levels.  Although  it  can  he  generalized  to  any  number  of  lewis.  ns 
complexity  increases  considerably  and  becomes  unleasable  in  multilevel  situations.  So. 
lor  one  and  two  dimensional  observations  which  have  more  than  two  intern::'.  ;c-. eS. 
we  need  to  search  lor  a  dillerent  approach. 

As  mentioned  in  the  first  section  of  this  chapter,  a  new  algoiithm  which  mes  the 
Kalman  1  iltcring  techniques  and  Dynamic  Programming  is  developed.  In  applications 
filtering  is  concerned  with  the  extraction  of  signals  from  noise.  If  the  observation  and 
original  scene  are  modeled  by  linear  stochastic  models,  a  solution  to  the  general 
filtering  problem  can  be  obtained  by  use  of  the  Kalman  Filter. 

Due  to  the  assumed  piecewise  characteristics  of  the  data,  the  realizations  of  Y 
can  be  modeled  by  the  state-space  models 


XfM  =  xt 


>t  =  Xt 


where  w(  is  Gaussian  zero  mean  i.i.d.  with  standard  deviation  a.  and  v  is  nonzero  at 
the  edges  between  regions  only.  By  defining  the  original  data  X  in  one  dimension  as 

X  =  {X(  .  t  €  LI,  L  =  (1.2 . N}.  the  one  dimensional  Markov  model  corresponding  to 

M  RF  has  joint  density 


lnP^(x)  =  P^g(Wi)  * 


and  the  likelihood  function  to  be  maximized  like  Equation  2.16  is 


fix)-  - -^-rV(yt-xt)2  +  pVg(Xi.Vl, 


where  the  first  term  on  the  right  hand  side  of  this  equation  comes  irom  the  noise  and 
the  second  term  is  the  potential  defined  before. 
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The  Kalman  Filter  tor  the  state-space  model  in  Equation  5.0  and  3.7  is  deiined 
by  the  recursion  Equation  3.1 1  as 

Vi =  \  ’  Kt-i^vv 

where  Kt  is  the  Kalman  gain.  When  an  edge  is  detected  the  gain  needs  to  be 
reinitialized  at  that  point.  To  do  this,  define  a  binary  sequence 


mt  =  vt 


tvhich  has  a  logic  "1"  at  the  edges  between  the  regions  and  a  logic  "0"  in  the  regions. 
Just  to  give  an  example,  let's  say  Tj  =  [tj.  t,  -  1]  and  T,  =  [t,  ,  t3  -  1]  where  T,.  T,  are 
adjacent  regions  and  t(  <  t,  <  t,.  We  can  define  V  j  =  vp  =  v  ,  =  1  as  shown  in 
Picture  T.l  In  ceneral 


1  if  t  -  l  g  Tj.  and  t  e  Tk  for  some  It 
0  otherwise 


V,9»1 


V,n  =1 


V,i  =1 


Figure  3.1  An  Example  for  Adjacent  Regions. 


Then  the  estimated  data  can  be  determined  from  the  noisy  sequence  v,  a»m:  E 
3. 10  and  3.1 1  as 


/V 

V 


K.-rii.YV 


where 


13. 


I 


and. 


r  1  m,  -  I 

T  ^  +■  1  otherwise 

Filtering  liquation  3.13  yields  the  estimate  of  x  as 


j=l 


t  i .  1 0 ) 


where  t  stands  for  the  distance  of  t  from  the  nearest  left  edge,  liquation  3. In 
expresses  the  Kalman  Filler  as  an  averaging  of  the  measurements. 

In  the  case  that  the  edges  are  not  known,  observe  the  following: 

1.  Given  the  noisy  data  v  for  any  binary  sequence  m,  the  estimate  \  is  well  defined 
by  Equation  3.13  -  3.15,  call  this  estimate  x(m); 

2.  The  g(.  ,  .)  term  in  the  Markov  model  Equation  3.S  and  the  likelihood  function 
Equation  3.9  can  be  interpreted  as  a  penalty  to  the  edges.  From  the 
dependence  of  x  on  the  sequence  nr  stated  above  it  can  be  written. 


g«xt 


1  if  m(  =  0 
- 1  if  mt  =  1 


1 3.1") 


Therefore  we  can  define  the  function 


Y<ny  = 


( .'.is i 


C  1  if  m(  =  0 
- 1  if  mt  =  l 

^  ^ 

For  any  j  e  L  let  y*.  xJ.  x1,  in1  be  the  sequences  and  y.  x.  x.  m  up  to  index  j.  for 

instance  v1  =  |vn,  v. .  v.|.  With  this  definition  the  likelihood  function 

liquation  3.9  associated  with  the  estimate  of  x  satisfies  the  recursion 


[i-^1(mj  +  1,=  ci(mi) 


v.  ,  ,-x. 
■  j +  i  i 


,(mJ 


*>r  +  Pv<mj  +  1) 


(3.iy» 


In  the  case  of  the  index  set  is  L  =  10,1 . X)  the  likelihood  function  liquation 

3.9  for  the  estimate  xim)  is 

flxtml)  =  (3.20) 

Bv  these  considerations  the  smoothing  problem  of  the  data  is  reduced  to 
maximization  of  the  likelihood  function  Equation  3.9  with  respect  to  all  binary- 
sequences  m.  This  can  be  done  using  the  recursive  formula  liquation  3.19  and 
Dynamic  Programming  techniques  given  the  details  in  Appendix  A. 

C.  SMOOTHING  OF  PIECEWISE  CONSTANT  SIGNALS  IN  TWO 

DIMENSIONS  USING  KALMAN  FILTERING 

There  is  a  problem  in  extending  Kalman  Filtering  techniques  from  one  dimension 
into  two  dimensions  due  to  the  lack  of  a  causal  state-space  model  for  higher 
dimensions.  Anyway,  we  can  determine  a  two  dimensional  recursive  formulation  like 
Equation  3.13  -  3.15  by  assuming  the  estimate  x  as  the  average  value  of  the  noisy  data 
within  the  regions. 

Proceeding  just  like  the  case  in  one  dimension,  a  smoothing  algorithm  can  be 
developed  to  maximize  the  likelihood  function  given  below: 

X  X  X 

Cc k.t)  =  V|y(k.t)-x(k.t)i:  +  P^g(x1,a.xka.1)  +  pVg(xk<l.xk.l  t)  (3.21) 

k.t  k.t  k.t 


To  do  this,  first  assume  the  edges  to  be  known,  and  define  M  =  !mk[!  where  ny  (  £ 
!e0.ej.e,.e,}  for  each  < k.t )  e  L.  e(  being  the  four  combinations  of  edges  corresponding 
to  the  four  possibilities  in  clique  ry  1  as 


These  combinations  of  edges  is  shown  in  Figure  3.2. 


Figure  3.2  Illustration  of  Edges. 

Some  definitions  are  given  below  referring  to  Figure  3.3  for  each  point  <Lti  6  L  to 
obtain  the  smoothing  algorithm: 

Definition  l :  rk  t  is  the  number  of  points  between  (k.t)  and  the’  closest  edge 
above. 

Definition  2:  zk  t  is  the  average  value  of  the  noisy  data  y  in  the  <  r,  .  x  1 1  region. 
Definition  3:  Xk  t  is  the  number  of  points  in  the  homogeneous  region  surrounded 
by  row  k.  column  t  and  the  line  of  edges. 

Definition  4:  xk  t  is  the  average  of  noisy  data  y  in  the  region  in  which  f  t  is. 
Using  the  definitions  given  above,  the  filtering  algorithm  can  be  developed  as 
following 


Zk.t  Zk-l,c  ~  K,t‘-k,t"Zk-l.i* 


k.t 


I  otherwise 


k 


k.t- 1 


Tk.t  tlnV,e 


tk  (  otncnvise 


1'his  algorithm  is  based  on  the  fact  tluu  tbr  any  disjoint  pair  of  regions  A  .  B  c  I 
average  ot  a  noisy  signal  denoted  by  v  on  the  regions  A  .  ft  .  A  U  B  are  relate.:  as 

m  A  U  B)  =  a |  \( A  i  a,\(B>  o'  ' 


with  a,  =  A  A  U  B’,  a,  =  B  A  U  B  where  .\  denotes  the  number  o:  eiemciU'  ,r 
the  set.  1  tie  recursion  liquation  3.24  computes  the  average  on  A  and  Equation  3.2< 
the  average  on  A  U  B.  As  in  the  one  dimensional  case  liquations  3.21  -  3.24  rep  re  sen1 

A 

a  well-delmed  mapping  which  associates  an  estimate  mMi  to  any  field  of  edges  \|  g 
;o().e1.e,,e,!  •  i  •  :  .  Analogous  considerations  as  lor  the  one  dimension..!  ..oc  !ea« 
to  the  likelihood  function  as 

t!\[mfvk[)]j  =  x  ■s'2[\|mivk  t)|'  '31< 


where 


r 


m,vk.t>  =  < 


M  ifvk,t  =  °o 
0lfvk.t  =  el  -e2 

l  -1  ifvk.t  =  e? 


So,  the  estimate  x  can  be  determined  by  searching  over  all  possible  sets  of  edges  M  i; 
order  to  find  x  which  maximizes  the  likelihood  I  unction  Equation  3.21.  Ile.vcve:  the 
Dvnamic  Programming  algorithm  is  used  to  ma\imi/e  the  likelihood  luiuuon  in 
dimensions. 


IV.  RESULTS  OF  SIMULATIONS 


A.  ONE  DIMENSIONAL  SIGNALS 

We  started  the  simulations  by  creating  a  one  dimensional  two  level  binary 
sequence  which  has  123  points  as  shown  in  Figure  4.1.  Gaussian  zero  mean  noise 
with  different  signal  to  noise  ratios  (SNRs)  calculated  with  the  formula  given  by 
Equation  3.2  was  generated  and  added  to  the  original  signal  by  using  1MSL  library 
subroutines.  The  noisy  sequences  are  given  in  Figure  4.2.  4.3,  4.4  and  4.5.  Here,  the 
standard  deviations  for  the  noise  are  25.  50.  "5  and  100  respectively. 


Figure  4.1  Original  Binary  Sequence. 


Figure  4.7  Segmented  Sequence  for  i;  =  So 


Figure  4.9  Segmented  Sequence  lor  <7=  1< hi  and  p  =  b.“2. 

For  the  case  of  data  with  unknown  levels  we  developed  a  different  segmentation 
algorithm.  The  difficulty  is  to  estimate  the  data  levels,  as  well  as  to  determine  the 
segmentation.  This  segmentation  algorithm  is  based  on  Kalman  Filtering  and 
Dynamic  Programming.  First,  this  method  has  been  applied  to  tne  one  dimensional 
noisy  signals  shown  in  Figures  4.2,  4.3.  4.4  and  4.3.  The  noisy  data  was  filtered 
twite  i  two  passes).  The  two  passes  results  are  shown  in  Figures  4. 1  i.i.  4.11  for  1=25 
and  3o. 


it 
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Figure  4.11  Segmented  Sequence  by  KF  with  <7=50. 

B.  TWO  DIMENSIONAL  IMAGES 

We  applied  the  two  dimensional  filtering  algorithm  to  the  test  images  which 
contain  three  different  regions  given  in  Figure  4.12.  In  this  test  image  there  are  four 
different  intensity  levels  (one  for  the  background  and  three  for  the  objects i.  Random 
Gaussian  noise  has  been  added  with  standard  deviation  <?  =  ID  and  <7=20  as  shown  in 
Figures  4.13,  4.14. 

The  result  of  filtering  is  shown  in  Figures  4.15  and  4. lb.  Notice  that  m  these 
cases  the  noise  has  been  removed  while  preserving  the  edges  between  the  regions.  The 
values  of  the  parameter  P  giving  the  best  results  are  P  =  2.o  for  <7=  10  and  P=  1.5  for 
<7=  20  showing  a  trend  of  this  algorithm.  The  value  of  p  depends  on  the  SNR  of  the 
noisy  picture.  In  particular  p  should  be  decreased  for  lusher  levels  of  noise. 


As  an  application  to  a  test  image  in  underwater  emironmcnt,  \vc  .sp-vivi  :d 
algorithm  to  a  512  x  512  picture  of  a  fish  corrupted  b\  additive  noise.  I  !.o  nay. 
image  is  shown  in  f  igure  4.17  and  the  filtered  one  is  shown  in  I  igu.e  4.1V  I  h.e 
significance  of  it  is  the  fact  that  after  filtering  'he  li<h  and  the  background  present  w  oil 
distinct  intensity  leads.  Although  applications  to  this  class  of  problems  are  still  under 
investigation,  the  fact  that  the  object  < the  fish  in  this  case)  presents  character. sties  well 
distinct  from  the  background  can  be  used  for  automatic  detection  or  recognition  m  an 


artificial  intelligence  framework. 


Figure  4.12  Original  l  est  Image. 


Figure  4.18  Segmented  Fish  with  <7=25  and  p  =  0.25. 


Analogously  we  tested  the  algorithm  on  a  16  level  checkboard  shown  m  Figure 
4.14  for  ditferent  levels  of  noise.  In  particular  the  values  of  the  16  levels  are  given  by 
100.  50.  ISO,  70.  200,  120,  60,  140,  75,  175,  90,  65,  130,  55,  190,  110  from  left  to  right 
for  each  line.  The  effects  of  additive  noise  are  given  in  Figures  4.20  and  4.21  for  a  -  10 
and  <7=20  respectively.  The  noise  is  always  assumed  to  be  Gaussian,  zero  mean  md 
white.  The  application  of  the  filtering  algorithm  is  shown  in  Figures  4.22  and  4.23 
using  [1  =  1.2  and  P  =  0.7  respectively.  As  expected  the  algorithm  filters  within  the 
regions  while  preserving  the  sharp  separation  between  adjacent  regions.  Also  shorn  n  in 
Figure  4.22  are  the  edges  between  regions,  as  detected  by  the  algorithm  which  shows 
the  reinitializations  of  the  Kalman  Filter  gains. 


V 


Figure  4.23  Segmented  Image  with  <j=20  and  p  =  0. 7. 
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V  CONCLUSION'S 


An  algorithm  to  segment  noisy  data  in  one  dimension  and  two  dimensions  has 
been  presented,  where  the  segments  are  piecewise  constant.  The  data  are  described  by 
state-space  models.  The  particular  feature  of  the  algorithm  is  to  search  lor  the  best 
sequence  of  edges  in  order  to  maximize  the  likelihood  function.  The  algorithm  has 
emphasized  piecewise  constant  data,  bur  it  can  be  extended  for  the  data  with  legions 
characterized  by  textures  to  which  we  associate  dillercnt  Autoregressive  iAR)  models. 
Typical  applications  are  not  only  image  segmentation  but  also  segmentation  in  speech 
analysis,  such  as  phenomenon  separation. 
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APPENDIX  A 

DYNAMIC  PROGRAMMING  ALGORITHM 


The  Dynamic  Programming  ns  an  algorithm  used  when  the  solution  to  a  problem 
may  be  considered  as  the  result  of  a  sequence  of  decisions.  An  optimal  sequence  of 
decisions  will  maximize  the  given  function,  in  the  case  that  is  used  in  this  stuJy,  the 
likelihood  function  fk. 

In  Dynamic  Programming  an  optimal  sequence  of  decisions  is  obtained  by 
making  explicit  appeal  to  the  Principle  of  Optimality.  In  the  cases  when  this  principle 
can  be  applied,  an  optimal  sequence  of  decisions  has  the  property  that  whatever  the 
initial  state  and  decision  are.  the  rest  of  the  decisions  must  make  up  an  optimal 
decision  sequence  according  to  the  state  resulting  from  the  first  decision.  Standard 
Dynamic  Programming  techniques  are  introduced  by  Bellman.  [Ref.  9| 

The  mathematical  equations  and  the  algorithm  is  given  below  for  the  case  that 
was  mentioned  in  the  second  section  of  Chapter  3.  By  defining  the  likelihood  function 
as 


\ 


lk'Vxi . V-Pl8<VYi>  ■  — I|>j-F(xi)|2 

j  =  0  2<T-  j  =  u 


or  in  recursive  form 


(A.  1 1 


+  1<VX1 . Xk+  fk<X0’Xl . ‘V  +  ^*Ts+  l'Xk+  l,Xk) 

where 


^k  +  l^xk+  l’V  P^xk+l’V'  ,  |Xk+ r^xk+ 1^ 


i  A.3 1 


and  setting  the  initial  condition  x^  =  0  and  also  assuming  that  xk  can  be  only  logic 
"<»"  or  "V.  we  can  determine  the  best  sequence  {xk  ,  k  =  0,....N)  which  maximizes 
f'qlXy . x>q).  For  this  purpose  define 
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yo)  =  max  fk(\() . xk  l.  u )  ,  v.4i 

<xu . Vi» 

and 

Vh  =  max  fk,xn . xk-i-  l>  ^ 

(  xij . \-i 1 

he  can  represent  the  problem  in  the  form  of  a  graph  as  in  figure  A.l.  with  X  -.tages 
where  the  nodes  at  the  k-th  stage  represent  the  state  x.  (either  "<f  or  'IX,  and  die 
branches  represent  the  gain  associated  to  the  transition  from  one  state  to  the  next.  In 
this  way  Jk(0)  and  1 )  represent  the  maximum  of  all  possible  gains  up  to  xk  =  q  and 
Xj.  =  1  respectively.  These  quantities  can  be  recursively  updated  as 

4+.  ,(0)=  max;Jk(0>  + ACk  +  ,(0.0)  .  Jk(  1  )+•  ACk  +  ,(l.0.)j  iA.6) 

and 

J,.  +  ,( 1 )  =  max !  Jk(0>  +  A£R  +  ,(0. 1  > ,  Jk(  1 )  +  ACk  +  ,i  1 , 1 )}  i  A. " > 

with  Afk  as  in  Equation  A. 3.  Also  we  can  keep  track  of  the  branches  which  yield  the 
maximum  likelihood  by  defining  Pointerk(0).  and  Pointer^  1)  at  each  stage  k.  1  he 

sequence  x  maximizing  ty.\0 . x^-)  is  therefore  obtained  by  backtracking  as  m  the 

following: 

Let  x^  be  such  that  Jyx^)  =  maxIJyO)  ,  jyi)J 
for  k  =  X  -  1  to  0 

\  =  Po,nterk  +  l<\  +  l) 
end  for 


APPENDIX  B 
COMPUTER  PROGRAMS 


i  ^  J 


m 


*  k  t  tt  k  r  *  t  -.  k  • 

T~ 

”  PROGRAM  B I MARY . DAT 

*  PURPOSE  To  generate  128  point  binary  sequence. 

*  REQUIRED  IMSL  ROUTINES  Hone. 

■k  ' 

*  IMPLEMENTED  BY  Lt.  Kani  HACIFASAOGLU  Jan.  .55? 

kkkkkkkkkkkk7tkk:xkkkkkk-xkki<-XAkk'K.-kkkk-k7t-k-kKkkkT:xr;7>;k'r7<;T;7Kr.  r.*r-r  r  »  ' 

INTEGER  IMITT , FIMITT ,F( 128 ) , TERM , SETBUF , K 
REAL  DWINDO , M07EA , DRAMA 

OPEN ( UNI T=1 , NAME= 1  BINARY . DAT  1 , ACCESS= ‘ DIRECT  1  , STATUS= ‘ NEW ' 
*  RECL=32 , MAaREC=123 ) 

DO  10  K=1 ,128 

IF (K  .GE.  0  .AND.  K  .LT.  30)  THEN 
F ( K ' =100 

ELSE  IF(K  .GT.  70  .AND.  K  .LT.  90)  THEN 
F  ( K  ,'  =  1 0  0 

ELSE  If f K  .GT.  120  .AND.  K  .LE.  128)  THEN 
F  f  K) =100 
ELSE 

F ( K ) =200 
END  IF 

WRITE ( 1 1 K )  (F(K) ) 

CONTINUE 
CLOSE (UNIT=1) 

CALL  INITT ( 430 ) 

CALL  TERM (2,1) 

CALL  SETBUF ( 2 ) 

CALL  DWINDO (-50.0,300.0,-50.0,400.0) 

CALL  MOvEA( -50 . 0 ,-50.0) 

CALL  DRAWA (300.0,-50.0) 

CALL  DRAWA (300. 0,400.0) 

CALL  DRAWA  -50.0,400.0) 

CALL  DRAWA ( -  50 . 0 , - 50 . 0 ) 

CALL  MOVEA  0.0,1000.0) 

CALL  DRAWA  0.0, 0.0) 

CALL  DRAWA  128.0,0.0) 

CALL  MOVEA(0. 0,0.0) 

DO  35  K=1 , 128 
K=K 
Y=F(K) 

CALL  DRAWA (X,Y) 

CONTINUE 

CALL  FINITT (0,0) 

STOP 

END 


+  fcjsxx^&  +  £7tJ*K,k"k’k-k'k'k’k’k'k'k-k'k -ft: -k  y:  jr -r.  yt-r.-r  +  irT: -TTt'k7'7t'i'k-*;7r  t: rrrrr  :Kr  -  t;-  - 

* 

*  PROGRAM  BINARY 2 . DAT 

:k 

*  PURPOSE  To  generate  one  dimensional  123  pom1;  bi 

*  sequence  '.;ith  zero  me  an  Gaus  ;i?.n"  noise  u. 

*  by  Dynamic  Programming  algorithm  to  segrv 

*  a  noisy  signal. 

X 

REQUIRED  IMSL  ROUTINES  GGNML 

*  IMPLEMENTED  BY  Lt.  Kara  HACIFASACGLU  Feb.  1987 

7C 

INTEGER  INI  IT  ,  FINITT , F( 123 >  .  TERM ,  SETB'JF  ,  I ,  MR  ,  L 
REAL  DWINDO ,  MOVEA ,  DF.AWA ,  GGNML  ,  R  ,  S  ( 123; 

DOUBLE  PRECISION  DSEED 
DSEED=65471 

OFEN’  ( UNIT-2 .NAME- 1 BINARY2 . DAT  1  , ACCESS* 1  DIRECT  1  , STATUS 
*  RECL=32 , MAXREC=123 ) 

DO  10  1=1,128 

IFfl  .GE.  0  .AMD.  I  .LT.  30)  THEN 
F ( I )=10Q 

ELSE  IF ( I  .GT.  70  .AND.  I  .LT.  90)  THEN 
F ( I ) =100 

ELSE  I F ( I  .GT.  120  .AND.  I  .LE.123)  THEN 
F ( I ) =100 
ELSE 

F ( I ) =200 
END  IF 
NR=  1 

CALL  GGNML (DSEED, NR, R) 

SIGMA=7  5 . 0 
P.=SIGHA*R 

R  s NOISE  FUNCTION  WHICH  WAS  GENERATED 
S ( I )=F( I ) tR 
WRITE ( 2 1 1 )  (S(I)) 

CONTINUE 
CLOSE (UNIT=2) 


CALL  IMITT (480 ) 

CALL  TERN (2,1) 

CALL  SETBUF  2) 

CALL  DWINDO (-50. 0,300. 0,-50. 0,400.0) 
CALL  MOVEA (-50.0, -50.0) 

CALL  DP.AWA  300.0,-50.0) 

CALL  DRAWA ( 300. 0,400.0 
CAIL  DRAWA ( -50.0,400.0) 

CALL  DRAWA (-50.0,-50.0) 

CALL  MOVEA(0. 0,1000.0) 

CALL  DRAWA ( 0 . 0 , 0 . 0 ) 

CALL  DRAWA  128.0,0.0) 

CALL  MOVEA (0.0, 0.0) 

DO  35  1=1,128 
X=I 

Y=S ( I ) 

CALL  DRAWA (X,Y) 

CONTINUE 

CALL  FINITT (0,0) 


vcvc 


FROGRAH  0MUR2 


PURPOSE  To  segment  3  noisy  binary  sequence  using 
Dynamic  Programming  algorithm. 


REQUIRED  IMSL  ROUTINES  Hone. 


IMPLEMENTED  BY  Lt.  Kani  HACIPASAOGLU  March  1967 


r-irx7*-^’Jc^x'^^,^7^K7^x^c^rx7T^:xx7,r+,xrrA^:x^:^cx-rrr7trcA^:^AAA^r^rc'r>t7r7ir-:K7,:>*r7crx;' 


INTEGER  N , K , FO , FI , OUT , PTR ( 128 , 0 : 1 ) , F ( 123 ) 

INTEGER  IN  ITT , FIN  ITT , TERM , 3ETEUF , J , M 

REAL  BETA , SIGMA , LNEWO , LUEW1 , MOVE A ,  LO  ,  LI ,  OUTS  ( 1 ;  12-3  )  ,  SUM  .  T -  E  ,  TOY 
REAL  DELVfl  :  123 , 0  : 1 ,  C  :  1 )  ,  DWINDO  ,  DRAMA ,  S  ( 123  )  ,  X ,  Y  ,  SSE  ,  SEE  ,  TEE 
S IGMA=1 50 . 0 
READ (  * ,  1 )  BETA 
FORMAT (F4. 1) 

FC=100 


OPEN ( UNIT=2 , NAHE= ' BIMARY2 . DAT  1 , ACCESS= 1  DIRECT ' , STATUS= ' OLD  1 
RECL=32 ,MAKREC=12B) 

DO  23  N=1 .123 
READ  (2'N  ( S ( M ) ) 

WRITE (  6  ,  * )  S(N),N 
CONTINUE 
CLOSE (UMIT=2) 

K=1 

L0=BETA- ( ( (S(K)-10Q)**2)/(2*(SIGMA+*2) ) ) 

Ll=-BETA-(( (S(K)-200)**2)/(2*(SIGHA**2))) 

DO  15  M=l,123 
K=M- 1 

DELV(K+1 , 0 , 0  )=BETA-  (((S(K+1)-F0)**2)/  (2*(  SIC-1IA**2 ) ) ) 

DELV (K+l ,1,0  =-BETA-  ( ( S (K+l ) -F0 ) **2 ) / ( 2*( SIGHA*"2 /  ) 

DELV  K+l,0,l)=-BETA-(  (S (K+l) -FI  ) **2 ) / ( 2*(SIGHA**2 ) } ) 

DELV (K+l, l,l)=3ETA-({  S (K  +  l ) -FI )  **2  )  /  ( 2*  ( SIGNAL ) ) ) 

IF ( (LO+DELV (K+l ,0,0))  . LT .  (Ll+DELV (K+l ,1,0)))  THEM 


( S (K) -200 ) **2 ) / 


DELV(K+1,0,0)=BETA-(((S(K+1)-F0)**2)/(2*(SIC-MA**2))) 
DELV (K+l ,1,0  =-BETA-  ( ( S (K+l ) -F0 ) **2 ) / ( 2+( SIGHA*"2 ;  ) 
DELV  K+1,0,1)=-BETA-  (S (K+l) -FI )  **2 ) / ( 2*(SIGHA**2 ) } ) 
DELV(K+l,l,l)=3ETA-((  S  (K  +  l ) -FI )  **2  )  /  ( 2*  ( SIGNAL ) ) ) 

IF ( (LO+DELV (K+l ,0,0))  . LT .  (Ll+DELV ( K+l , 1 , 0 )) )  THEM 
LMEW0=L1+DELV(K+1 ,1,0) 

PTR (M , 0 )=1 
ELSE 

LNEW0=L0+DELV(K+1 ,0,0) 

PTR(M, 0)=0 
END  IF 

IF((L0+DELV(K+1,0,1))  .LT.  (L1+DELV(K+1 , 1 , 1 ) ) )  THEN 
LNEW1=L1+DELV ( K+l , 1 , 1 ) 

FTR(M, 1 )=1 
ELSE 

LNEW1=L0+DELV (K+l ,0,1) 

PTR(H, 1 )=0 
END  IF 

L0=LNEW0 

L1=LNEW1 

CONTINUE 

TRACE  BACK  PROCEDURE 

OPEN (UNIT=3,NAME=‘ OUTS. DAT1  , ACCESS= 1  DIRECT  1 ,STATUS='HEW 
*  RECL=32 ,MAXREC=128) 

IF ( LO  .LT.  LI)  THEN 
J=1 
ELSE 
J=0 
END  IF 
OUT= J 

DO  44  K=127 ,0,-1 
IF (OUT  .EO.  0)  THEN 
OUTS (K+l )=100 
ELSE 

CUTS (K+l )=200 
END  IF 


OUT=PTR( K+l , OUT) 

WRITE (6,*)  OUTS (K+l ) , K+l 
CONTINUE 
DO  55  K=1  128 
WRITE(3'K) (OUTS (K ) ) 
CONTINUE 
CLOSE ( UNIT=3 ) 


CALL  INITT (480 ) 

CALL  TERIK  2,1) 

CALL  SET3UF ( 2 ) 

CALL  DWINDO (-50.0,300.0,-50.0,400.0) 
CALL  I10VEA  -50.0,-50.0) 

CALL  DRAWA  300.0,-50.0) 

CALL  DRAWA  300.0,400.0) 

CALL  DRAWA (-50. 0,400.0) 

CALL  DRAWA ( -50.0 , -50.0) 

CALL  HOVEA  0.0,1000.0) 

CALL  DRAWA  0.0, 0.0) 

CALL  DRAWA (150.0,0.0) 

CALL  HOVEA (0.0, 0.0) 

DO  3  K=1 , 128 
X=K 

Y=OUTS (K) 

CALL  DRAWA (X,Y) 

CONTINUE 
CALL  FINITT (0,0) 


A 


Sw 


PROGRAM  TEST.DAT 


PURPOSE  To  qenerate  a  test  image  which  has 
three  different  intensity  levels. 


REQUIRED  IMSL  ROUTIIIES  Hone. 

IMPLEMENTED  BY  Lt.  Kani  HACIPASAOGLU  May  1937 


BYTE  A( 123 , 123) 

INTEGER  R,XC1 , YC1 ,XC2 , YC2 , FI ( 1 : 128,1 : 123) ,F2(1 : 123,1 : 123) 


R=20 

XC1=60 

YC1=60 

YC2=60 

XC2=90 


OPEN (UMIT=5 , NAME- 1  TEST . DAT 1 , ACCESS= 1  DIRECT 1 , STATU5= 
RECL=32 (HAXREC=128) 


DO  10  1=1,128 
DO  20  J=1 , 123 


FI (I , J)=( ( (I-XC1)**2+(J-YC1)**2)-(R**2) 
F2(l/J)  =  (((I-XC2)*'i:2+(J-YC2)*,'2)-(R:t*2) 

( ( FI ( I , J )  .LE.  0)  .AND.  ( F2 (I , J )  .GT.  0) 


IF( ( FI ( I , J )  .LE.  0)  .AND.  (F2(I,J)  .GT.  0))  THEN 
A ( I , J  )  =  50 

ELSE  IF ( ( FI ( I , J )  .GT.  0)  .AND.  (F2(I,J)  .LE.  0))  THEN 
A ( I , J ) =1 50 

ELSE  IF ( ( FI ( I , J )  .LE.  0)  .AND.  (F2(I,J)  .LE.  0))  THEN 


A ( I , J ) =1 50 

ELSE  IF ( ( FI ( I , J )  .LE.  0)  .A 
A( I , J )=100 
ELSE 

ACl , J)=200 
END  IF 
CONTINUE 

WRITE ( 5 1  I )  (A(I , J) ,  J=1 , 128) 
CONTINUE 
CLOSE (UNIT=5) 

STOP 

END 


SNKRKKKf 


U 


Li<unjh»Li  M.tf*  I.*  m  j|«  _**».♦»*  j  i  i  , 


★  *******£x*x********:*x*:k*x*;fc**;fr;fr***:frx***:*r7‘c:4:***7i:7r7i'r7»-7r;'i':x:r;*r:- 

★ 

*  PROGRAM  TESTIM.DAT 

X 

*  PURFOSE  To  generate  a  test  image  with  zero  mean 

*  Gaussian  noise  to  be  segmented  by 

*  Kalman  Filter. 

* 

*  REQUIRED  IMSL  ROUTINES  GGNML 

*  IMPLEMENTED  BY  Lt.  Kani  HACIPASAOGLU  May  1937 

x 

x  A  *  *  x  *  x  x  x  *  x  *  *  x  x  *  *  x  x  *  *  x  *  *  x  *  x  *  x  x  A  *  *  *  *  *  x  x  x  *  *  x  x  x  x  x  x  x  x  x  x  x  x  x  x  * 

BYTE  A( 123 , 123) 

CHARACTEP*50  TESTIM 

INTEGER  MCI , YC1 , XC2 , YC2 , FI ( 128 , 128) 

INTEGER  NR, Rr, SIGMA, CR, F2( 123, 123) ,AA( 123,123) 

REAL  GGNML, R 

DOUBLE  PRECISION  DSEED 

DSEED=65471 

WRITE (*,100) 

FORMAT ( 'ENTER  IMAGE  DATA  NAME') 

READ (*, 101) 

FORMAT ;A50) 

CR=20 

KC1=60 

YC1=60 

YC2=60 

XC2=90 

OPEM(UNIT=l , NAME= ' TESTIM.DAT ' ,ACCESS= ' DIRECT ' ,STATUS=' 
*  RECL=32 , MAXREC=128 ) 

WRITE (*,222) 

FORMAT ('ENTER  SIGMA1) 

READ (*,333)  SIGMA 
FORMAT (12) 

DO  10  1-1,128 
DO  20  J=1 , 128 

FI ( I , J)=( ( ( I-XC1 )**2+( J-YC1 )**2) - (CR**2 ) ) 

F2(I , J)=( ( (I-XC2)**2+( J~YC2)**2)-(CR**2) ) 

IF ( ( FI ( I , J )  .LE.  0)  .AND.  (F2(I,J)  .GT.  0))  THEN 
A( I , J )=50 

ELSE  IF ( ( FI ( I , J )  -GT.  0)  .AND.  (F2(I,J)  .LE.  0))  THEN 
A( I  ,  J)  =  150 

ELSE  IF(( FI ( I , J)  .LE.  0)  .AND.  (F2(I,J)  .LE.  0))  THEN 
A ( I , J ) =100 
ELSE 

A ( I , J )=200 
END  IF 
NR-1 

CALL  GGNML (DSEED. NR, R) 

Rr=SIGMA* jnint (R) 

AA( I , J )=A( I , J)+Rr 
IF (AA(I , J )  .LT.  -128)  THEN 
AA( I , J )=AA( I , J ) +256 
ELSE  IF ( AA( I , J)  .GT.  127)  THEN 
AA( I , J)=AA( I , J) -256 
ELSE 

AA ( I , J ) =AA ( I , J ) 

END  IF 

A ( I , J ) =AA ( I , J ) 

CONTINUE 

WPITE(l’I)  (A(I,J),J= 1,128) 

CONTINUE 
CLOSE (UNIT=1) 

STOP 

END 


k  — 

*  FROGRAM  FINAL . DAT 

*  - 

*  PURPOSE  To  generate  a  test  image  winch  lias  16 

*  different  regions.  * 

k  r 

*  REQUIRED  IMSL  ROUTINES  None.  * 

"T  ^ 

*  IMPLEMENTED  Ex  Lt.  Kani  HACIPASAOGLU  June  1957 

kkkkkkkkkkkkkkkkkkkkkTrkkkkkk  kkkkT^^kkrcTr-xkkkkkkk  xkTKkkkkk-r-K^ck  r.-xTC 


BYTE  A (123, 123) 

OFEN(UNIT=l , NANE= 1 -INAL . DAT  1  , ACCE3S= 1  DIRECT  1  , STATUS* 
■  RECL=32,UAUEEC=123; 

DO  10  K=1 . 32 
DO  20  L=1 , 32 

A  ( K ,  L )  =  1 0  0 
A  ( K , L+32 ) =50 
A(  K , L+64 ;  =  130 
A(K,L+96)=70 

A (K+32 , L )=200 
A(K+32,L+32)=120 
A  Kt32,L+64)=60 
A i K+32 , L+96 )  =  140 

A( K+64 , L )=75 
A(K+64,L+32)=175 
A<K+64(L+64)=90 
A(K+64 , L+96 ) =65 

A(K+96 ,L)=130 
A(K+96 ,L+32)=55 
A(K+96,L+64)=190 
A(K+96 ,L+96 )=110 
CONTINUE 
CONTINUE 
DO  30  1=1,128 

WRITE ( 1 1  I )  (A( I , J) , J=1 , 123) 

CONTINUE 
CLOSE (UNIT=1) 

STOP 

END 


I 

a 


R 

i 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 

101 


222 

333 


C 

C 


20 

10 


90 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk'kkkkkk-kk-xk^kk  kkk-k  xir  r 

k  r 

*  PROGRAM  NFINAL.DAT 

k  ■*■ 

PURPOSE  To  generate  a  test  image  which  has  10  ~ 

*  different  regions  with  zero  mean  Gaussian 

*  noise  to  be  segmented  by  Kalman  Filter.  * 

*  REQUIRED  IHSL  ROUTIilES  GGMML  * 

k 

*  IMPLEMENTED  BY  Lt.  Kani  HACIFASAOGLU  June  1937 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk  kkkkkkk  kkkk-r:kr. 


BYTE  A( 123, 123) 
CHARACTER* 50  MFINAL 


INTEGER  HR, Rr, SIGMA, AA( 128, 128) 
REAL  GGMML , R 
DOUBLE  PRECISION  DSEED 
DSEED=65471 
WRITE (*, 100) 

FORMAT  1  ENTER  IMAGE  DATA  NAME’ 
READ (",101) 

FORMAT (A50) 


) 


OPEN ( UMIT=1 ,NAME=' MFINAL3.DAT' , ACCESS= 1  DIRECT  1 , STATUS* 1  HEW 1 , 
*  RECL=32 ,MAKREC=128) 

WRITE (*,222) 

FORMAT/  1  ENTER  SIGMA'  '> 
kEAD  v  * , 333  )  SIGMA 
FORMAT ( 12 ) 

DO  10  K=l,32 
DO  20  L=l,32 

A ( K , L ) =100 
A  K,L+32)=50 
A(  K , L+64 )  =  180 
A(K , L+96 )=7Q 

A(K+32 ,L) -200 
A(K+32,L+32)=120 
A(K+32,L+64)=60 
A(K+32 , L+96 )=140 

A ( K+64 , L)=7  5 
A(K+64,L+32)=175 
A(K+64,L+64)=90 
A ( K+64 , L+96 ) =6 5 

A(K+96,L)=130 
A  K+96,L+32)=55 
A(K+96, L+64  =190 
A( K+96 , L+96 ) =110 
CONTINUE 
CONTINUE 
MR=1 

DO  80  1=1,128 

DO  90  J=1 , 128 

CALL  GGMML (DSEED, NR, R) 

Rr=SIGMA*inint(R) 

AA( I , J)=A( I , J ) +Rr 
IF(AA(I,J)  .LT.  -123)  THEN 
AA( I , J )=AA (I , J ) +256 
ELSE  I F(AA ( I , J )  .GT.  127)  THEN 
AA( I , J )=AA( I , J ) -  256 
ELSE 

AA ( I , J ) =AA ( I , J ) 

END  IF 

A ( I , J ) =AA ( I , J ) 

CONTINUE 

WRITE ( 1 1  I )  (A(I , J) , J=1 , 123) 
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*  PROGRAM  CRISTI2D  * 

*  K 

'■  PURPOSE  To  segment  a  two  dimensional  noisy  image 
using ’Markov  Random  Field,  Dynamic 
'  Programming  and  Kalman  Filtering  techniques.  ' 

*  T 

*  REQUIRED  IHSL  ROUTINES  None.  * 

7T  - 

*  IMPLEMENTED  BY  Lt.  Kani  HACIPASAOGLU  Aug.  1987 

common  ny,  v ( 512 , 512 ) , tau( 512 ) , xtop ( 512 ) , iedge , ymax 
integer  start,  dir,  cycle 


character*2Q  image,  imout 
byte  s ( 512 , 512 ) ,  v(512,512) 
write!*,  777 ) 

777  format  ('  Enter  Input  Image  File  Name-.') 

read  (*,  778)  image 
write!*,  782) 

732  format('  Enter  Output  Image  File  Name:') 

read!*,  773)  imout 
773  format(a20) 

open  (2,  name=image,  access= 1  direct 1  ,  status= ' old1  ) 
c  get  data  from  file 

write  (*,  779) 

779  format( 1  Enter  Image  Size1) 

read!*,  780)  npoints 
730  format(i3) 

write (*,  5551 

read  !*,  666)  sigma,  beta 

•.•rite  (*,  790) 

790  format('  see  edges?  (yes=l)') 

read  !*,  791)  iedge 

791  format(il) 

write!*,  810) 
write (*,  792) 

792  format!'  Enter  scale  factor  (1.0=no  scale;1) 

310  formats'  output  yout=yth+scale*(y-y th)  1  ) 

read!*,  793)  scale 

793  format ( f 10 . 2 ) 

write!*, 811) 

811  format('  1:  yth=(yav+ymax)/2 ,-  2 : yth=0 . 9ymax ;  3:  enter  ytiv;  ??') 
read!*,  791)  iyth 
if  (iyth.eq.3)  then 
write(*,  794) 

794  format!'  Enter  yth:‘) 

read  (*,  793)  yth 
endif 
c 

do  5  k=l,  npoints 

read(2'k)  (S(k,j),  j=l,  npoints) 
do  6  j=l,  npoints 
temp=s(k, i ) 
if  ( temp . It . 0 . 0)  then 
temp=temp+256 
endif 

y(k, j )=temp 
6  continue 

5  continue 

close  (unit=2) 
c 
c 

555  format!'  ENTER:  sigma,  beta') 

666  format(2fl0 .4) 

c 
c 


***  main  loop  *** 


667 


wr lte ( * ,  667) 

format ( 1  Start  Computinq  ...') 
ymax=Q.O 
yav=0 .0 
ncount=0 
nx=npoints 
ny=npoints 

gam=l  .0/(2. 0*sigma'A:*2) 

***  initialize 

do  30  j=l,ny 
tau( j )=0.0 
xtoo  ( j )=0 . 0 
continue 

nx=nx-l 
***  main  loop 
do  40  i=l,nx 
start=l 
dir=+l 
cycle=l 

call  line(i,  start,  dir,  cycle,  gam,  beta,  rov;av) 

start=ny 

dir=-l 

cycle=2 

call  line(i, start,  dir,  cycle,  gam,  beta,rowav) 
continue 

***  compute  average  of  output  image 
yav=(ncount*yav+rowav)/ (ncount+1 ; 
ncount=ncount+l 

write(*,  668) 

format^'  •••  End  Computing.  Mow  transfering  data  to  disk 
***  output  data  *** 

if  (iyth.eq. 2)  then 
yth=ymax*0.9 
endif 

if  ( iyth.eq. 1)  then 
y th= ( yav+ymax ) / 2 . 0 
endif 

write(*,  800) 

formate  yave,  ymax,  yth  =') 
write (*,  801)  yav,  ymax,  yth 
format(3(x,  fl0.4)) 


do  60  k=l ,  npoints 
do  61  j=l,  npoints 

ys=yth+scale*(y(k, j )-yth) 

if  (ys .gt. 255.0)  then 
ys=255 .0 
endif 

if  (ys.le.0.0)  then 
ys=0 .0 
endif 

itemp=jnint(ys) 
if  (itemp.gt . 127)  then 
itemp=itemp-256 
endif 

v(k, j )=itemp 
continue 
continue 


c 


50 

c 


90 

c 

c 


999 

c 

c 

c 

c 


110 

c 


open  (unit=3,  name=imout,  access= 1  direct 1 
status= 1  new 1 ,  recl=128,  maxrec=512) 


do  50  j=l  npoints 
write(3 1 j )  (v(j,k),  k=l, npoints) 
continue 
close  (unit=3) 


stop 

end 


subroutine  line(i,  start,  dir,  cycle, 
irt,  air, 


gam , 


beta, rowav) 
xtop(512 ) iedge  ,  ymax 


integer  start,  dir,  cycle 
common  ny,  y(512,512),  tau(512) 
real  xhat(4,512),  lambda(.4,512) 
real  ttau(4) ,  txtop(4),  tlambda(4),  txhat(4),  te(4) 
real  newtau (4 , 512 ) ,  newxtop(4 , 512) 
real  e(4, 512) 
integer  pointer (4 , 512) 


***  scan  from  START  by  DIR  (=+l,  -1) 
j=start 


***  initialize  first  column 
do  90  k=l,4 

xhat  ( k , start )=y( i , start) 
lambda (k, start) =0.0 
e ( k, s tart )=0 . 0 
continue 


rowav=0 . 0 
ncount=0 
i=j+dir 

***  state=^left  edge,  upper  edge) 


0=no  edge,  l=edge 


***  state  (0,0) 
do  110  k=l .  4 
ttau(k)=tau( j )+l 

txtcp( k)=xtop( j )+(l . 0/ ttau(k) )*(y ( i, j )-xtop( j ) ) 
tlambda(k)=lambdaik, j-dir )+ttau(k) 
z=txtop(k)-xhat(k, j-dir) 
txhat(k)=xhat(k, j -dir)+( ttau(k)/tlambda(k) )*z 


:l=(y(i, j+dir ) - txhat(k) )**2 
:2=(y(i(i )-txnat(k) )**2 


e  1  =  < 

e3=( i+1 ’ j )-txhat(k|)**2 


te(k)=e(k, j-dir)Tgam*(el+e2+e3)/3  -  2*beta 


if  (k.eq.l)  then 

emin=te(k) 

endif 

if  (te(k) .le.emin)  then 
im=k 

emin=te(k) 

endif 

continue 


newtau (4 , j )=txtop( im) 
lambda (4 , j )=tlambda( im) 
pointer (4 , j )=im 
newxtop(4, j )=txtop(im) 
xhat (4, ] )=txhat(im) 
e (4 , j )=te(im) 


***  state  (0,1)  *** 
do  115  k= 1 , 4 

tlambda(k)=lambda(k,  j-dir)+l 
z=y(i, j )-xhat(k, j-dir) 
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txhat  (k)=xhat(k,  j  -  dir )  +•(  1 . 0/t lambda (k)  )*z 

el=(y(  i,  j+dir)- txhat  (  k)  )**2 
e2=(y (i , 3 )  -  txhat ( k) > **2 
e3=(y(i+l , j )-txhat(k) ) '*2 
te(k)=e(k, j -dir )+gamx(el+e2+e3)/  3 

if  (k.eq.l)  then 

emin=te(k) 

endif 

if  ( te(k) . le . emin)  then 
im=k 

emin=te ( k) 
endif 

continue 

newtau( 1 ,  j )  =1 
ne«xtop( 1 ,  -j  )=y(  i ,  j  ) 
lambda ( 1 , j )=tlambda( im) 
pointer ( 1 , j )=im 
xhat ( 1 , 3 )= txhat ( im) 
e  ( 1 ,  j )  =  te ( im) 


***  state  (1,0) 
do  120  1;=1 ,4 
ttau(k)=tau( j )+l 

txtop(k)=xtop(j)  +  ( 1 . 0/ ttau(k) )*(y (i, j ) -xtcp( j ) ) 

el=(y(i, j+dir )-txtop(k) )**2 
e2={V[i,j )  -txtop(  it)  )**2 
e3=(y( l+l , i )-txtop(k) )^*2 
te(k)=e(k,  j-dir)  +  gamrt(el+e2+e3)/3 

if  (k.eq.l)  then 

emin=te(k) 

endif 

if  ( te(k) . le . emin)  then 
im=k 

emin=te(k) 

endif 

continue 


newtau(2 , j )=ttau( im) 
nev;xtop(2(]  )=txtcp(im) 
lambda(2 , j ;=tau( j ) 
pointer ( 2 , j )=im 
xhat (2, ] )=txtop(im) 
e (2 , j )=te ( im) 


***  state  (1,1) 
newtau(3 , j )=1 
newxtop ( 3 , 1 ) =y ( i , j ) 
lambda(3 . j )=1 
xhat(3, j)=y(i,j) 

e3=(y(i+l,3)-y(i,j)  )K*2 
z2=gam*(el+e3)/3 

if  (e(4,  j-dir) . It. e( 1 , j-dir ) )  then 
e ( 3 , 3 )=z2+e ( 4 , j -dir )+beta 
pointer(3, j )=4 

else 


e(3,j)=z2+e(i , j-dir)+beta 
pointer ( 3 , j )=1 

endif 


if  ( j .gt . (1-dir ) . and. 3 . It . (ny-dir ) )  then 
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goto  999 
endif 

***  backtrack  *** 

emin=e(l , j ) 
do  150  k=l , 4 

if  (emin.le.e(k, j ) )  then 
emin=e (k, j ) 
ini=k 
endif 


150 

c 

continue 

383 

c 

y(i, j)=xhat(im, j) 

c 

write (*,  111)  im 

111 

format(i3) 

if  (cycle. eq. 2)  then 

xtop(j  )=nev;xtOD(  im  .  j) 
tau( j )=newtau^im , ] ) 

***  compute  max  intensity 
if  (y(i, j ) .at.ymax)  then 
ymax=y  ("i , ] ) 

endif 


**  row  average 

rowav=(ncount*rowav+y(i , j ) )/ (ncount+1 ) 
nccunt=ncounttl 

mark  the  edge  with  a  black  entry 

if  ( im. ne .4. and. iedge .eq. 1 )  then 
y(i,j)=o.o 
endif 
endif 

im=pointer(im, j ) 
j=j -dir 

if  ( j .gt. (1-dir) .and. j .It. (ny-dir) )  then 
goto  388 

endif 

return 

end 


m 
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